#include <math.h>

#ifndef M_PI
#define M_PI 3.1415926535897932385
#endif

double besselj0 (const double x) {
    double t, t2, a0, p0, q0, ta0, bj0;
    double pi = M_PI;

    if (x == 0.0) {
        bj0 = 1.0;
    }
    else if ( (0.0 < x) && (x <= 4.0) ) {
        t=x/4.0;
        t2=t*t;
        bj0=((((((-.5014415e-3*t2+.76771853e-2)*t2-.0709253492)*t2+.4443584263)*t2 -1.7777560599)*t2+3.9999973021)*t2-3.9999998721)*t2+1.0;
    }
    else {
        t=4.0/x;
        t2=t*t;
        a0=sqrt(2.0/(pi*x));
        p0=((((-.9285e-5*t2+.43506e-4)*t2-.122226e-3)*t2+.434725e-3)*t2-.4394275e-2)*t2+.999999997;
        q0=t*(((((.8099e-5*t2-.35614e-4)*t2+.85844e-4)*t2-.218024e-3)*t2+.1144106e-2)*t2-.031249995);
        ta0=x-.25*pi;
        bj0=a0*(p0*cos(ta0)-q0*sin(ta0));
    }
    
    return bj0;
}

